Operational analysis for COVID-19 testing: Determining the risk from asymptomatic infections

Testing remains a key tool for managing health care and making health policy during the coronavirus pandemic, and it will probably be important in future pandemics. Because of false negative and false positive tests, the observed fraction of positive tests—the surface positivity—is generally different from the fraction of infected individuals (the incidence rate of the disease). In this paper a previous method for translating surface positivity to a point estimate for incidence rate, then to an appropriate range of values for the incidence rate consistent with the model and data (the test range), and finally to the risk (the probability of including one infected individual) associated with groups of different sizes is illustrated. The method is then extended to include asymptomatic infections. To do so, the process of testing is modeled using both analysis and Monte Carlo simulation. Doing so shows that it is possible to determine point estimates for the fraction of infected and symptomatic individuals, the fraction of uninfected and symptomatic individuals, and the ratio of infected asymptomatic individuals to infected symptomatic individuals. Inclusion of symptom status generalizes the test range from an interval to a region in the plane determined by the incidence rate and the ratio of asymptomatic to symptomatic infections; likelihood methods can be used to determine the contour of the rest region. Points on this contour can be used to compute the risk (defined as the probability of including one asymptomatic infected individual) in groups of different sizes. These results have operational implications that include: positivity rate is not incidence rate; symptom status at testing can provide valuable information about asymptomatic infections; collecting information on time since putative virus exposure at testing is valuable for determining point estimates and test ranges; risk is a graded (rather than binary) function of group size; and because the information provided by testing becomes more accurate with more tests but at a decreasing rate, it is possible to over-test fixed spatial regions. The paper concludes with limitations of the method and directions for future work.


Introduction
Entering the third year of the 2019 coronavirus disease (henceforth COVID-19) pandemic, it is clear that the world was woefully underprepared, in many different ways, for dealing with it. ItThe current pandemic is an illustration of the natural evolutionary play [1] so that one should expect future pandemics and lose no time in preparing for them while dealing with the present one.
Some authors have suggested that it is appropriate to prepare for the next pandemic as one prepares for war [2][3][4][5]. Operational analysis grew out of the scientific approach to operational questions in World War II [6][7][8][9]. One of the key tenets of operational analysis is to model the process as well as the data [7,10]. Among the advantages, when the process is modeled, one knows the true state of world, which allows assessment of the quality of the analyses by comparison of analytical outcomes with a known situation. This gives confidence that the methods will work when the true state of the world is unknown.
Models for dynamics and control of the disease, prioritizing hospital care, and setting policy [11][12][13][14][15][16][17][18] require information about the health status of the population. This is determined by testing for infection, which thus emerged as a crucial component of managing health policy during the current pandemic and will probably be key in future pandemics [19]. For example, the Global Influenza Surveillance and Response Network (the "flu network" [20,21]) established in 1952 played a key role in the early responses to the COVID-19 pandemic. The time is now to prepare for future testing.
Testing is complicated an individual is in an early stage of the infection may give a false negative test, infected individuals may be asymptomatic and thus not tested, and symptomatic but uninfected individuals may give false positive test results. Thus, test errors involve both false negative tests, and false positive tests, in which an uninfected individual tests positive [22][23][24][25][26][27]. These are called the false negative probability [28], denoted here by p FN , and false positive probability, denoted here p FP . It may one day be possible to drive the false positive probability to zero with improved specificity of tests, but the ontogeny of the disease within an individual means that there will always be false negative tests [25,26].
A starting point for the interpretation of testing results is to envision that a population is divided into infected (antigen positive) and uninfected (antigen negative) individuals, with the goal of estimating the fraction of infected individuals (the incidence rate) from the number of positive results P when T tests are given. Because of both kinds of test errors, the surface positivity rate P/T (which is observed; henceforth simply called positivity rate) will generally differ from the incidence rate (which is not observed). It is natural and intuitive to ask for the unobserved incidence rate that is most likely given the test results; this is called the Maximum Likelihood Estimate (MLE) of the incidence rate.
Brown and Mangel [29] and Mangel and Brown [30] (also see [31,32] where similar methods are used) show that the Maximum Likelihood Estimtate for the incidence rate, denoted bŷ which is to be interpreted asf ¼ 0 if the right side of Eq 1 is negative. As will be explained below application of Bayesian methods allows determination of a probability distribution for the incidence rate when the right side of Eq 1 is negative. In addition to a point estimate for the incidence rate, it is valuable to have a range of incidence rates that are consistent with the model and the data since then one can bound the incidence rate and its associated risk of further infection (Eqs 3 and 4 below). That is, forecasting for a pandemic can be improved by using a predictive distribution, rather than the point estimates in Eq 1, (cf. [33]).
Rangeðf Þ is symmetrically distributed around the true range with very small mean error between the two [30], so that lower and upper limits for the estimated incidence rate arê f lower ¼f À 0:5 � Rangeðf Þ andf upper ¼f þ 0:5 � Rangeðf Þ.
Mangel and Brown [30] also show how likelihood methods can be used to obtain a test range when positivity is 0 (so thatf ¼ 0). In this paper, we will show how to determine the test range when 0 < P/T < p FP However, for the remainder of this section, we assume that P/T > p FP so that the estimate of incidence rate is strictly positive. Eqs 1 and 2 lead to the operational recommendation that one should stratify testing data according to test errors. When this is not possible, one should stratify tests according to the estimated time since exposure, assign best estimates to the test errors, and conduct sensitivity analyses of the results.
Test results can play an important role in policymaking because they can be used to determine the risk of spreading infection associated with groups of different sizes. Doing so requires a definition of risk. We define the risk to be including at least one infected individual in a group of specified size. The risk Rðh;f Þ associated with a group of size h when the estimate for incidence rate isf is [29,30] Eq 3 allows us to explore the risk ramifications of groups of different sizes. Were the true incidence rate known, we replacef by f t (Fig 1) Replacing the estimate of incidence rate by its maximum and minimum values f max ¼f þ 0:5 � Rangeðf ) and f min ¼f À 0:5 � Rangeðf Þ in Eq 4 allows us to bound the acceptable group size consistent with the level of acceptable risk.
In Fig 2, I show 16 realizations of acceptable group size using the simulation methods described in [30]. The dotted line shows the group size consistent with the level of acceptable risk when the incidence rate is f t , the solid black line is the group size using the estimate in Eq 4, and the red and blue lines are the group sizes using the maximum and minimum estimates for incidence rate, f upper ¼f þ 0:5 � Rangeðf ) and f lower ¼f À 0:5 � Rangeðf Þ, respectively. One key observation is that the group size determined if the incidence rate were known (dotted line) falls between those determined from the upper and lower limits of incidence rate determined by the test range.
The first purpose of this paper is to show how to obtain the test range when there is no information on symptoms and positivity is less than the probability of a false positive test. The second purpose of this paper is to generalize Eqs 1-4 and develop the analogue of When there is information on symptoms (Fig 3), a fraction f t of the population is infected and symptomatic; a fraction ρ t f t is infected and asymptomatic; a fraction g t is uninfected but symptomatic; and the remaining fraction, 1 − f t (1 + ρ t ) − g t , is neither infected nor symptomatic. Infected individuals have probabilities of a false negative test, denoted by p SFN and p AFN where the subscript S and A correspond to symptomatic and asymptomatic individuals, respectively. Uninfected individuals have probabilities of a false positive test denoted by p SFP and p AFP , respectively. Using testing information, we seek point estimates and the analogue of test ranges for the unknown incidence rates and ratio of asymptomatic to symptomatic cases.

Determining test range with no information on symptoms and positivity less than the probability of a false positive test
In this case, as with Eqs 1-4, there is a single unknown incidence rate, which we continue to denote by f t . The methods used are generalized when there is information on symptoms, so this section is a warm-up to the harder problem.
The probability of obtaining a positive test when the incidence rate is f is The first term on the right hand side of Eq 5 corresponds to individuals who are infected and have a true positive test; the second term corresponds to individuals who are not infected and have a false positive test. Sixteen realizations of the group size as a function of acceptable risk using the simulation methods described in [30]. In all panels, the number of tests is T = 2500, the true incidence rate is f t = 0.05, and the probabilities of false negative and false positives tests are 0.25 and 0.05. The dotted line shows the group size consistent with the level of acceptable risk when the incidence rate is f t , the solid black line is the group size using the estimate in Eq 1 determined using the positivity rate from the individual realization of the simulation, and the red and blue lines are the group sizes using the maximum and minimum estimates for incidence rate, f upper ¼f þ 0:5 � Rangeðf ) and f lower ¼f À 0:5 � Rangeðf Þ, respectively. One key observation is that the group size determined if the incidence rate were known (dotted line) falls between those determined from the upper and lower limits of incidence rate determined by the test range. https://doi.org/10.1371/journal.pone.0281710.g002 When T tests are given, the number of positive tests P is binomially distributed with parameters T and p + (f) [30][31][32], which we write as P ¼ Bð�; T; p þ ðf t ÞÞ, where BðP; T; p þ ðf ÞÞ ¼ T P À � p þ ðf Þ P ð1 À p þ ðf ÞÞ TÀ P . The likelihood of an incidence rate f given the test data P and T has the same form [30][31][32], but is a function of the incidence rate f conditioned on the values of the test data In S1 Section in S1 File, we show that the maximum likelihood estimate for the fraction of the population infectedf satisfies Solving this equation forf gives Eq 1. When P � Tp + (f), we setf = 0 and arrive at the nettlesome case of this subsection.
In   We convert from a likelihood to a probability distribution by assuming a uniform prior on f and use Bayes's theorem to write the probability density for f given the test data (also see S2 Section in S1 File): Although the denominator in Eq 8 can be written in terms of the classical beta function [48], it is most simply viewed as a constant obtained by using a very fine discretization of the interval [0, 1].
When the maximum of the likelihood occurs at the boundary f = 0, the probability φ(f) will also have its maximum at the boundary. In this case, the test range is no longer symmetrical but is an interval [0, f 0.95 ], where f 0.95 is the value of incidence rate such that R f 0:95 f 0 ¼0 �ðf 0 jP; TÞdf 0 ¼ 0:95 (or the equivalent when a summation instead of an integral is used in Eq 8).

Analysis when there is information on symptoms
The operational situation with information on symptoms. We assume that T tests are administered to a population in which some individuals are symptomatic and others are not (recorded at the time of testing) and each individual tested has either a positive or negative test for coronavirus. As described in Fig 3, there are now four classes of individuals: • A fraction f t of the population is symptomatic and infected (antigen positive); these individuals have a probability of a false negative test p SFN .
• A fraction g t of the population is symptomatic but not infected; these individuals have a probability of a false positive test p SFP .
• A fraction ρ t f t of the population is infected but not symptomatic; these individuals have a probability of a false negative test p AFN .
• The remaining fraction of the population, is neither infected nor symptomatic; these individuals have a probability of a false positive test p AFP .
When this situation holds, three kinds of test data are generated: • The number P of positive tests.
• The number T S of symptomatic individuals.
• The number P S of symptomatic individuals who tested positive.

Point estimates for the fractions of infected and uninfected symptomatic individuals and the ratio of asymptomatic to symptomatic infected individuals.
The following causal chain characterizes the operation of testing: • The total number of tests, T, leads to number of symptomatic individuals in the sample, T S .
• T S leads to the number of positive tests of symptomatic individuals, P S .
• T, T S , and P S combined lead to the remaining positive results, P − P S of T − T S tests from asymptomatic individuals.
As above, we let Bð�jN; pÞ denote a binomial distribution with number of samples N and probability of a positive event p, where the dot can run from 0 (no positive event) to N (only positive events). If p S denotes the probability of sampling a symptomatic individual, p S+ the probability of obtaining a positive test from a symptomatic individual, and p A+ the probability of obtaining a positive test from an asymptomatic individual, the test results have distributions The probabilities on the right sides of Eqs 9-11 are constructed from the assumptions summarized in Fig 3. The fraction of individuals who are symptomatic is f t + g t so that Since p S+ is the probability that an individual tests positive given that the individual is symptomatic, from the definition of conditional probability The probability that an individual tests positive given that the individual is asymptomatic is computed similarly: We letf ;ĝ andr denote the maximum likelihood estimates (MLEs) for the fraction of the population that is infected and symptomatic or not infected and symptomatic respectively, and for the ratio of the fraction that is infected and asymptomatic to that which is infected and symptomatic.
When a random variable has the binomial distribution Bð�jN; pÞ, given K positive events, the MLE for p isp ¼ K=N (see S1 Section in S1 File) so that the MLEs for the probabilities in Eqs 9-11 are T S /T, P S /T S , and (P − P S )/(T − T S ) from which we concludê frð1 À p AFN Þ þ ð1 Àf ð1 þrÞ Àĝ Þp AFP ¼ ð1 Àf Àĝ Þ P À P S T À T S : Eqs 15 and 16 are independent ofr and can be rewritten aŝ which we write asf ¼ c 1 ðP S ; T S Þĝ , where c 1 (P S , T S ) is the combination of terms multiplyingĝ on the right side of Eq 18 and the test errors are suppressed. With this notation, we substitute into Eq 15 and solve to obtainĝ Thus, bothf andĝ are known; they are random variables because P S and T S are random variables.
We now rewrite Eq 17 aŝ frð1 À p AFN À p AFP Þ ¼ ð1 Àf Àĝ Þ P À P S T À T S À p AFP ð1 Àf Àĝ Þ; ð20Þ let c 2 = 1 − p AFN − p AFP , and solve forr to obtain Since the right side of Eq 21 depends on the test data T S , P S , and P,r is also a random variable.
Eqs 18, 19 and 21 generalize Eq 2 to the case in which symptomatic and asymptomatic individuals are identified at the time of testing. We have already thus generalized the method in [29,30] to obtain point estimates of the fractions of the population of infected and symptomatic, uninfected and symptomatic, and infected and asymptomatic individuals. We next explore the properties of these point estimates and then generalize the notion of test range and compute and the risk of groups of different sizes including asymptomatic infected individuals.
The means of the estimates. We compute the means of the estimates, continuing to use f t , g t , and ρ t to denote their true values, with two goals. We explore 1) whether the estimates in Eqs 18, 19 and 21 are unbiased, in the sense that their expectations (over the stochastic sampling process) are the underlying true values generating the data, and 2) if there is a bias how to characterize it.
The mean ofĝ . We begin by rewriting Eq 19 aŝ In S2 Section in S1 File, I show that 1 þ cðP S ; T S Þ ¼ 1À p SFN À p SFP 1À p SFN À ðP S =T S Þ so that we can rewrite Eq 22 asĝ Since the denominator in Eq 23 is a constant, the expectation ofĝ is In the S2 Section in S1 File, I show that EðT S Þ ¼ Tðf t þ g t Þ and EðP S Þ ¼ T½f t ð1 À p SFN Þ þ g t p SFP �, from which it follows that Eðĝ Þ ¼ g t ; the expected value ofĝ is the true value that underlies the testing process.
The mean off . We begin by multiplying the top and bottom of the right side of Eq 18 by T S to obtainf and now use the version ofĝ on the far right side of Eq 23 to obtain Taking expectations on the far right side of Eq 26, we obtain so that the expected value off is the true value that underlies the testing process. The mean ofr. We begin with Eq 21, rewritten aŝ Eq 28 is a nonlinear function off andĝ and involves the quotients of the random variables. We can approximate the expectation ofr using the delta-method [30,49], which involves Taylor expansion of the right hand side of Eq 28 to second order and then taking expectations. (Details are in the S2 Section in S1 File). The result is where Var(X) and Cov(X, Y) denote the variance and covariance of random variables X and Y, which arise from the second order Taylor expansion. The right side of Eq 29 shows that the leading term in the expected value ofr is the true value generating the data and that this is corrected by variances and covariances that account for the nonlinearity in Eq 24.

Joint properties off andr via likelihood analysis
Eq 2 for the test range can be obtained by direct manipulation of the relevant random variables [30]. When we separate symptomatic and asymptomatic infections, the compatibility interval for the incidence rate is replaced by a compatibility region (CR) for the incidence rate of symptomatic infected individuals and the ratio of asymptomatic to symptomatic individuals. Because Eqs 18 and 19 are nonlinear in the test results (which are random variables) and Eq 21 is also nonlinear inf , the analytical approach used in Mangel and Brown [30] is less feasible now.
We develop the analogue of Eq 2 for test range by using likelihood analysis [50][51][52], exploiting the general property that for a smooth and well-behaved likelihood (which those that follow are), a 95% CR can be approximated by finding the range of variables for which the loglikelihood is below the peak log-likelihood by 1.96 times the number of free parameters. This is essentially a generalization of the Gaussian approximation to the binomial distribution [53] that leads to Eq 2 [30].
We denote the test results byT S ;P S , andP. For any values of f, g, and ρ, the rules of conditional probability imply (suppressing the dependence on T which is known) Each term on the right side of Eq 30 is a binomial distribution. In particular, for any values of f, g, and ρ, ð31Þ the probabilities on the right side of Eqs 31-33 are, respectively, p S , p S+ , and p A+ in Eqs 12-14 for any values of f, g and ρ, rather than their true but unknown values. When data T S , P S , and P are obtained, the likelihoods, given the data, that the state of the environment is f, g, and ρ are The likelihoods in Eqs 37 and 38 are products of binomial distributions that are well approximated, for sufficient numbers of tests, by the appropriate Gaussian distribution [30,53]. In the results, we will explore log-likelihoods for both the binomial distributions and their Gaussian approximations.
Simplifying the likelihoods. Keeping our eyes on the prize of computing the risk of including infected but asymptomatic individuals in groups of different sizes, we focus on f and ρ when constructing the CR. Exploring the likelihood is more convenient if one can eliminate having to deal with g explicitly. Two methods are the profile likelihood and the marginal likelihood [49]; both reduce the number of parameters from 3 to 2.
For the profile likelihood, we replace g in Eqs 37 and 38 by the MLEĝ , so that By numerical exploration, I found that for the operational questions modeled here, the two methods give virtually the same results for the answers. Were we interested in the tails of the likelihood, this might not be the case. Since the profile likelihood is computationally much speedier, I report results using it. The third Rscript in S4 Section in S1 File allows one to explore the differences between marginal and profile likelihoods for the symptomatic data.
The compatibility regions from the profile likelihood. I computed the approximate 95% CR from the total profile likelihood using a generalization of the method of Hudson [51] by first finding the maximum value of the profile log-likehood and then determining the region in f, g, or f, ρ-space in which the log-likelihood was 2 � 1.96 = 3.92 below its maximum value.
I did computations using R Studio 1.0.143 with underlying R 3.6.1 GUI 1.70 El Capitan build (7684) on an iMac running Mac OS 12.1.

Test range with no information on symptoms and positivity less than the probability of a false positive test
When positivity is less than the probability of false positive test, φ(f|P, T) has, similar to the likelihood, its maximum at f = 0 and monotonically declines. The peak value of φ(f|P, T) and the rate of decline depend on the positivity and the number of tests (Fig 5).
These normalized likelihoods In Fig 5 have a test range that depends on the number of tests (Fig 6). As with the situation in which positivity exceeds the probability of a false positive test, the test range declines with test numbers but at a decreasing rate.

Point estimates, compatibility regions, and risk when there is information on symptoms
In the base case for Monte Carlo simulations, I set N = 1000 replicates of T = 1500 tests. Since simulation and test errors scale as the reciprocal of their values, these choices have inherent errors of the order of 3%, which are sufficient to understand the qualitative patterns and most of the quantitative patterns. I chose the parameters for the true state of the world and the test errors from those reported in [22][23][24][25][26][27]: f t = 0.05, g t = 0.04, and ρ t = 1.5 and the test errors are p SFN = 0.25, p SFP = 0.03, p AFN = 0.5, and p AFP = 0.003. S3 Section in S1 File contains results for other choices of the true but unknown state of the world.
For the likelihood calculations and the associated risk computations, I first assume that the test results are the expected values T S , P S , and P, which is a reasonable assumption when T is large enough, after which I allow the test results to vary more widely. For the base case parameters, the mean values are � T S ¼ 135; � P S ¼ 58:05, and � P ¼ 118:0575. Since actual test data can only produce integer values, I rounded the � P S and � P to 58 and 118, respectively. Doing so gives the point estimatesf ¼ 0:04995;ĝ ¼ 0:04005, andr ¼ 1:50119 (significant digits included to illustrate how little accuracy is lost by the rounding process; since the true values are f t = 0.05, g t = 0.04, and ρ t = 1.5).
Illustrative simulated data. The n th replicate of the simulation of the testing process yields estimatesf n ;ĝ n , andr n . In Fig 7, I show the first 100 values of the simulation replicates.
Each circle represents the value off n ;ĝ n , orr n on the n th replicate of the simulation. The thick red lines represent the averages over the entire 1000 simulations. There are also black lines at the true values of the three parameters.
The means off n andĝ n essentially sit on top of the true values, as we would expect from the analysis in Eqs 22-27 showing that Eðf Þ ¼ f t and Eðĝ Þ ¼ g t . To quantify this agreement, I computed the mean relative error (ME) for the three estimates. For example, forf , it is For the simulation illustrated in Fig 7, MEðf n Þ ¼ 0:0039 and MEðĝ n Þ ¼ 0:0036 (i.e., both a fraction of a percent).
The lower right panel of Fig 7 has an expanded y-axis to show that the mean of ther n exceeds ρ t . For this run of the simulation, MEðr n Þ ¼ 0:0141. While this is less than 2.0%, it is almost four times larger than the mean errors off n andĝ n .
In Fig 7, the thin dotted lines show the means of the estimates ±1.96 times their standard deviations. These are a naive 95% compatibility interval under a Gaussian approximation because they ignore the other two parameters. Even so, for the full set of 1000 replicates of the simulation, the fractions of points outside this naive interval are 0.045, 0.055, and 0.05, respectively, forf n ,ĝ n , andr n .
We conclude that the formulas for the MLEs accurately capture their true values. The specific results, of course, depend on the simulation results and the number of tests given (both addressed in the next section). For example, in a different run of 1000 simulations of the testing process, the mean relative errors were -0.0073, -0.0006, and -0.033 forf n ,ĝ n , andr n respectively, and the fractions of points outside of the naive 95% CR were 0.054, 0.035, and 0.045 for f n ,ĝ n , andr n , respectively.
Likelihood, compatibility regions, and the risk of groups of different sizes. In order to focus on a single value of "test data" we continue using the expected values of T S , P S , and P. After exploring the situation when test data are the mean value, we will vary the test data.
The likelihood of the symptomatic data. On the way to the goal of estimating the fraction of asymptomatic infections, it is worthwhile to briefly stop and explore the likelihood of the symptomatic data, which are independent of ρ (Eq 37). In Fig 8, I show the likelihood when the means of T S and P S are the test results. In this figure, the white dot denotes the true values of parameters and sits at the peak of the heat map.
When the incidence rate is f and the fraction of symptomatic individuals who are uninfected is g, the mean value of T S = T(f + g) so that we expect a negative correlation between values of f and g, which is evidenced in the figure by the orientation of the contours of likelihood. For the purposes of the risk calculation, the most important role of the likelihood function of the symptomatic data is to provide the MLE value of g for construction of the profile likelihood in Eq 40, to which we now turn.
The likelihood of all the data. In Fig 9, I show the profile likelihood (Eq 40) for f and ρ when the test data are the mean values of of T S , P S , and P. The banana shape of the contours of the 95% CR computed is a result of the nonlinearity in Eq 28. In this case, the likelihood is centered at the true values of the parameters (shown by the white dot), and when Eq 28 is converted to a function ρ(f) by using the MLEĝ and replacingf by f andr by ρ, the true values of the parameters sit on the resulting curve, which runs through the middle of the 95% CR.
One property of the range formula in Eq 2 is that test range declines as 1= ffi ffi ffiffi T p . That is, although the range declines as the number of tests increases, it does so at a decreasing rate [30,]. This observation is more than an academic point, because it has the operational implication that it is possible to over-sample by providing too many tests in a single spatial region (also [30]).
In Fig 10, I explore the consequences of simultaneously increasing the number of tests and relaxing the assumption that the test data are the mean values of T S , P S , and P so that the true values (the white dots) no longer sit in the middle of the 95% CR or on the curve ρ(f) and the contours move in space, as determined by the test results. As with Eq 2, contours shrink as the number of tests increases, but at a decreasing rate. The prize: The risk of including asymptomatic infected individuals in groups of different sizes. We are now able to compute the risk of including asymptomatic but infected individuals in groups of different sizes, to generate a curve analogous to Fig 1. For any values of f and ρ, the fraction of asymptomatic infected individuals in the population is ρf. Hence, the analogue of Eq 3 is In Fig 11,

Discussion
It is important to recognize that the analysis presented in this paper is a procedure that allows one to go from test information to the risk of including infected individuals in groups of various sizes when there is no information on symptoms at the time of testing or to the risk of including asymptomatic individuals in groups of different sizes when there is information on testing. Rather than being binary (risky or not), this risk is graded and the specific details of the relationship between group size and risk depends upon the operational details of testing such as test numbers and errors. Once these are specified, the procedures can be employed. Let us now consider three limitations of the methods developed here. First note that Eq 26 has the same problem as Eq 1 when the positivity is very small. To see this, we factor out T S on the far right side of Eq 26 to obtainf so that if the positivity rate among symptomatic individuals falls below the probability of a false positive test among symptomatic individuals,f is less than zero. As in the situation with no information on symptoms, the operational interpretation is that we then setf to 0. Alternatively, we may generalize the analysis for the simpler case by putting a prior on the parameters and determining the CR in that manner. Second, an objection may be made that the binomial distribution underlying the analysis relies on the strong assumption that tests are independent events but that often groups of people will test together so that modeling the testing process requires an aggregated distribution. This is a fair objection, however: 1) the binomial distribution is the appropriate starting point, and if the sample is large and diverse enough (e.g., from many different testing sites), the independence assumption should be at least approximately valid; and 2) a negative binomial distribution of the form used in ecology to model aggregated counts [48, pp. 103-111] is a natural starting point for extending the work here.
Third, an objection may be made that we have assumed the values for test errors rather than estimating them. DiCiccio et al. [54] show that estimating test errors at the same time as incidence rate is a much more complex problem, and is likely one whose solution is not easily transferred to recommendations for practice. An alternative is to stratify test results by both symptomatic or not and time since putative exposure, have approximate values for the test errors for each time since exposure, and conduct sensitivity analysis by varying the values of the test errors. Furthermore, Mangel and Brown [30, pp. 25-27] show how to generalize Eq 1 for the case of a distribution of test errors using the delta-method. A similar extension of Eqs 18, 19 and 21 is a potential next step in this work.
In some locations, individuals are already asked whether they are symptomatic or not at the time of testing. For example, Nomi Health in Utah requires a self-reporting form for obtaining a coronavirus test, and the form includes yes or no questions such as: "Do you have a fever, a cough, new or increased shortness of breath, decreased smell or taste, a sore throat, muscle aches or pains, a headache, congestion or a runny nose, nauseas or vomiting, diarrhea, fatigue?" During the 2020-2022 academic years, natural experiments in testing were occurring on college campuses [55]. The results of those tests will provide a trove of information to explore with the methods developed here.

Conclusions
In conclusion, the approach of modeling and simulating the process of testing before analyzing testing data leads to a range of insights and at least the following operational recommendations: • At the time of testing, collect information on whether and individual is symptomatic or not.
• At the time of testing, collect information on putative time since exposure to infection.
• Conduct experiments to obtain information on means and variances of test errors.
There is much to be done and no time to lose before the next pandemic.
Supporting information S1 File. Including a brief review of the binomial distribution and likelihood, a mathematical appendix with details of calculations in the main text, sensitivity analysis when there is information on symptoms, and codes that generate the results in the main text and sensitivity analysis. (PDF)